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ABSTRACT 

We present an algorithm, AROFAC2, which detects the (CP- 
)rank of a degree 3 tensor and calculates its factorization into 
rank-one components. We provide generative conditions for 
the algorithm to work and demonstrate on both synthetic and 
real world data that AROFAC2 is a potentially outperform- 
ing alternative to the gold standard PARAFAC over which it 
has the advantages that it can intrinsically detect the true rank, 
avoids spurious components, and is stable with respect to out- 
liers and non-Gaussian noise. 

Index Terms — Tensor Decomposition, Tensor Factoriza- 
tion, Approximate Algebra, Simultaneous Diagonalization 

1. INTRODUCTION 

Polyadic decomposition of tensors into their canonical com- 
ponents (= canonic polyadic resp. CP-decomposition) and de- 
termining the number of those (= the rank) is a multidimen- 
sional generalization of the Singular Value Decomposition 
and the matrix rank, and a reoccurring task in all practical sci- 
ences, appearing many times under different names; first dis- 
covered by Hitchcock [ 1 1 and then re-discovered under names 
such as PARAFAC or CANDECOMP 0, it has been ap- 
plied in many fields such as chemometrics, psychometrics, 
and signal processing HE] ID- An extensive survey of many 
applications can be found in J7JE1- 

Considerable effort has been devoted to develop theory 
and methodology for the CP-decomposition, however many 
fundamental issues are still unresolved. The mathematical 
theory concerning CP-decompositions of tensors which are 
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not matrices is only partly understood; also, while there exist 
several methods to calculate the CP -decomposition of a ten- 
sor J9l [TO), they are extrinsical in the sense that a structure- 
agnostic loss function is optimized and also highly sensitive 
to outliers or non-Gaussian noise - problems which have been 
heuristically attempted to cope with (e.g. ifTTI '). Moreover, 
determining the rank of a noisy tensor remains a problematic 
task despite the existence of heuristics lfl2l . 

In this paper, we present AROFAC2, a method for cal- 
culating the CP-decomposition of a low-rank degree 3 ten- 
sor and its rank, which is based on theoretical considerations 
and intrinsical calculations making use of the algebraic struc- 
ture of degree 3 tensors, part of which have already surfaced 
in |fl3l , Specifically, we show how the algebraic structure can 
be used to obtain components one-by-one by alternating pro- 
jections - a technique which draws inspiration from lfl4ll - and 
how to reduce determination of rank to a clustering problem. 
Due to its structure-awareness, our algorithm only finds the 
numerically stable components while avoiding spurious ones, 
and determines the correct rank; it is also less sensitive to 
outliers or noise. We demonstrate its superiority to existing 
approaches on synthetic data and a chemometrics data set. 

2. TENSOR FACTORIZATION AND 
SIMULTANEOUS SVD 

We briefly review the basic definitions of rank-one tensor de- 
composition, and introduce some notation. 

Notations 2.1. Tlie set of {n\ x n 2 x n^-tensors of degree 3 
is denoted by C" lX " 2Xn3 . For A e C" 1 x ™ 2 x ™ 3 , the matrices 

Ax, ... , A„ 3 with A k = (aijk)i<i<m 

l<j<n 2 

are called the 3-slices of A. 

Definition 2.2. Let A e C" lXn2X " 3 . Then, a decomposition 

r 

A = ^Ui®Vi®Wi with u t e C" 1 , Vi £ C™ 2 , w 4 G C" 3 , 

i=l 



is called a rank r canonic polyadic decomposition (or CP- 
decomposition) of A. The Ui,Vi,Wi are called (rank-one- 
)components of A. The U{ are called mode-1-, the Vi are 
called mode-%, and the w^ mode-3-components. The tensors 
Ui ® Vi C3> Wi are called {rank-one-)f actors. 

The (CP-)rank of A, denoted by rk(A), is the smallest r 
such that A has a rank r CP-decomposition. 

This paper proposes an algorithmic solution for the fol- 
lowing problem: 

Problem 2.3. Let A £ C" lX " 2Xn3 . Determine r = rk(A) 
and a rank r CP-decomposition of A. 

and for its approximate version, i.e., the case where A is 
noisy but of low tensor rank, and one wants to find a CP- 
decomposition of the noiseless A. 

Related to that is the problem of finding a simultaneous 
singular value decomposition (SVD): 



Assumption 3.1. A £ C niX " 2X " 3 ,rk(A) < min(m,n2,n 3 ). 

It is probably known or folklore that in this case the CP- 
decomposition is unique, nevertheless we were not able to 
retrieve an exact reference. Thus, we provide a proof instead 
(for the notion of genericity, consult the appendix of 1 15 1): 



n± X 712 X ": J , 



Theorem 3.2. Let A £ 

rk(A) = r < min(rii, ri2, n^), let 



be generic with rank 



A : 



i=l 



u, <g> Vi <g> Wi with m G C™ 1 , Vi G C" 2 , Wi G C" 3 



be a CP-decomposition of A. Then any CP-decomposition of 

A can be obtained by replacing the Ui,Vi, Wi by 

XiUi, PiVi, A^~ vT Wi, where Xi, Vi G C x for 1 < i < r. 



Proof. Keep the notation from Proposition 2.5 including A 



Problem 2.4. Let A 1} ... , A n , G 



and Mj . Due to the observation in the proof of Proposition 2.5 



r, matrices U £ 
Wi,...,W n3 G 

1 < k < n 3 . 



",V G 



i xn 2 Determine a rank 
and diagonal matrices 



<r such that A k = U ■ W k ■ V 1 for all 



that (i) and (iii) differ only in notation, the statement of this 
theorem is equivalent to proving that there is only one unique 
way to present the Ai as 



In fact, Problems 2.3 and 2.4 are equivalent. We briefly 
prove that and add one additional characterization which will 
be important in the application: 



for all 1 < i < 7J 3 . 



Proposition 2.5. Let A G 



ini x 7i2 x ri3 



let 



A, 



i ^4-n 3 G 



with rank one matrices Mj and numbers Ay, up to renumber- 
ing and the obvious rescaling by replacing Mj with p,jMj and 



be the 3-slices of A. Then, the fol- \j with pJ 1 M J . Now since A is generic, the Mj are generic 

rank-one matrices, and the A<y are generic numbers. Thus, in- 
terpreting the Ai as rows of a (77,3 x n\ n-z) -matrix A, the Mj 
as rows of a (r x niri2)-matrix M, and the \j as elements 
of a (n 3 x r)-matrix A. The presentation above can be re- 
formulated as A — KM. The assumption < r and Propo- 
sition 
Ai 



lowing are equivalent: 

(i) A is a tensor of rank r in the sense of Definition ^.. 2\ 

(ii) The Ai have a simultaneous SVD of rank r in the sense of 
ProblemW4\ 

(iii) There are rank one matrices M%, . . . , M r G C" lXrt2 and 
coefficients Xij G C, 1 < i < n.3, 1 < j < r, such that 

r 

Ai = XijMj for all 1 < i < 113. 

Proof. (i)4=>(ii): There is only a notational difference which 
results from writing the columns of U as Uj, the columns of 
V as Vi, and, for i fixed and k running, the vectors formed by 
the i-th diagonal entries of the W% as Wi, and the A k as the 
3-slices of A. 

(i)<S=Kiii): Again, the difference is only notational. Write 
Mj = Ui^Vi and Ay the i-th component of Wj. □ 

3. UNIQUENESS OF CP-DECOMPOSITION AND 
LINEAR FACTORING 

By the definitions above, it is not a-priori clear whether the 
CP-decomposition, including terminology such as component 
or factor, is well-defined, in particular whether it is unique. 
The method we present in this paper exploits certain proper- 
ties of tensors with low rank. The main assumption on the 
tensor A to factor is the following: 



2.5 (iii) thus imply that the Mj lie in the span of the 
Now since the Ay are generic, A is a completely generic 
matrix. Thus, a different presentation of A would correspond 
to the existence of an invertible (r x r) matrix P such that 
A = kP~ 1 M' with a matrix M' whose rows correspond to 
rank-one-matrices, i.e., M' = PM. But a linear combination 

r 

Mi = Y / PijM j 

i=l 

of the (ni x 712) matrices Mj has rank one if and only if 
exactly one of the py (with i fixed) is non-zero, since r < 
min(ni, 712). Since P is of full rank, this implies that P is 
the product of a (r x r) permutation matrix with a full rank 
diagonal (rxr) matrix. But this is, as stated above, equivalent 
to the statement to prove. □ 

Thus, under our assumptions components are unique up 
to scaling and numbering, and factors are unique up to num- 
bering. 



Proposition 2.5 together with the uniqueness guarantee in 



Theorem 3.2 gives the following statement: 
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Proposition 3.3. Let A £ c niX " 2X " 3 be generic with rank 
rk(A) = r < n 3 , let A u ...,A n3 £ C" lX ™ 2 be the 3- 
slices of A. Then, up to scaling, there are r unique vectors 
\V>, AW g C™ 3 such that 

"3 

M(AW)=g\ ( * ) A j 

i=l 

/zas ra«£ one. Moreover, /or eac/z fc, let M(A^) = Uk eg) V). 
be the ( rank- one -)SVD. Then, the Uk and Vk are mode-1- and 
mode-2- components of A, belonging to the same factor. 

4. THE AROFAC2 ALGORITHM 

We propose an algorithm which computes rank and CP- 
decomposition of a degree 3 tensor A. It uses the routine 
FindRankOne which finds mode-1- and -2-components 
which we will present first as Algorithm [T] In step [T] the 
tensor is first decomposed into 3-slices Aj. In step [2] an ap- 
proximate representation V for their span is calculated. This 
can be a PCA of the A, i.e., principal values or components, 
or a numerical span of lower dimension (if the true rank of A 
is known, the dimension should equal the rank). In step [3] a 
matrix M E V is randomized. This can be a random matrix 
in an exact span, or a matrix which is, e.g., sampled from a 
Gaussian with covariance matrix follows the estimated sam- 
ple distribution of the A. Then steps [5] and [6] are repeated 
until convergence of M is attained. Step [5] takes (a possibly 
non-square) M to its third power, magnifying its largest sin- 
gular value and diminishing the others. Step [6] projects M 
onto V and normalizes the result. Projection can be achieved 
by exact projection, or re-scaling, e.g., according to the prin- 
cipal values of the Aj- After convergence is reached, M will 
be approximately of rank one, and its first singular vectors, 
which are obtained in step [8] will be estimates for a mode-l- 
and a mode-2-component of A. 

Algorithm 1 FindRankOne(A) Input: A £ C" lXn2X ™ 3 
Output: One random mo de-\- component u € C" 1 and one 
random mode-2-component v £ C™ 2 of A in the same factor. 

1: LetAi, .... 

2: Calculate 
span(Ai, . . 

3: Randomize M £ V. 

4: repeat 

5: M<-M- M T ■ M 

6: M<-V V (M) 

7: until M has converged 

8: Calculate approximate rank-one SVD of M = u ■ v T 
9: Return u, v 

Algorithm [2] uses Algorithm [T] to obtain a full CP- 
decomposition and an estimate for the rank of A. In step [T] 
several candidate estimates for components of all modes are 



obtained. Mode-3-components can be obtained by switching 
coordinates in A (e.g., switch the first with the third). If the 
mode-3-components are not relevant for the problem at hand, 



e.g., in the setting of simultaneous SVD as in Problem 2.4 



A n3 be the 3-slices of A. 
an approximate representation 

, Ar,„ ) . 



V of 



this can be omitted. In steps [2] and step [3] a clustering al- 
gorithm is applied to estimate the number of cluster centers 
and cluster the candidate components. This can be done by 
different methods, or one single algorithm. Our implemen- 
tation uses the mean shift algorithm which also estimates 
the number of cluster centers 11611 . Since Algorithm [T] links 
pairs of components, this information can then be applied in 
step|4]to the estimated cluster centers in order to link pairs to 
full triples (which is unnecessary if only the first two modes 
are considered), e.g., by majority vote or vote weighted by 
closeness of a pair to the cluster center. The corresponding 
decomposition is then presented as the estimated solution in 
step [5] 

Algorithm 2 AR0FAC2 (A) Input: A £ C" 1 x ™ 2 x ™ 3 
Output: Rank and approximate CP -decomposition of A. 
1: Repeat FindRankOne(A) to find sets 51,52,53 of po- 
tential mode-1-, -2-, and -3-components of A. 
2: Use clustering algorithm to determine number r of clus- 
ters for Si , S2 , and £3 
3: Cluster Si,S2, and S3 to obtain cluster centers 
U\, . . . ,u r of S\, centers v\, . . . , V r of 52 and Wi, . . . , w r 
of S 3 . 

4: Use the information from FindRankOne(A) to renum- 
ber the Ui , Vj , Wk such that for each I, the components 
ui,ve,we belong to the same factor. 

5: Return r as the rank of A and A = Y2i=i u i ® v i ® Wi 
as its CP-decomposition. 



5. EXPERIMENTS 



First, we demonstrate our algorithm on simulated toy data. 
The input tensor consists of 3-slices A& = ^21=1 ^ik u i V J 



compare Proposition 2.5 Each slice is generated as follows: 
Exact singular vectors Ui £ K." 1 , w, £ W 12 , 1 < i < r are 
sampled independently and uniformly from the rt!-sphere and 
n 2 -sphere, respectively. The A,^ are sampled independently 
and uniformly from the standard normal distribution. Then, 
to each matrix Afc, 1 < k < n^, noise is added in the form 
of a (ni x n2) matrix whose entries are independently sam- 
pled from a normal distribution with mean and covariance 



e £ M + . Figure 1(a) shows the accuracy of the estimated U 
and V for n\ — 50, ri2 = 60, 713 = 70, r = 10, e = 0.1 for 
PARAFAC (top row) vs AROFAC2 (bottom row). The esti- 
mation quality is essentially the same, however for AROFAC2 
the correct rank r = 10 has been detected automatically. 

In a second numerical evaluation, we analyze the noise 
robustness of the AROFAC2 algorithm. The data is generated 
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(a) estimated coefficients (b) estimated rank 



Fig. 1. (a) Absolute values of the correlation coefficients 
between the true and estimated components for PARAFAC 
(top row) vs AROFAC2 (bottom row). The correct solution 
has the same support as a permutation matrix, (b) Estimated 
rank for increasing noise levels e. True rank is r = 10. 

as before, but the noise level e is increased from e — 0.01 to 
e = 0.6. In Figure [T(b)| we see, that the correct rank r = 10 
has been found for e < 0.35 and is mainly overestimated for 
larger noise levels. 

Finally we apply our algorithm to a publicly available 
data set from chemometrics. The Dorrit fluorescence data 
ifTTl fTBl contains 27 synthetic samples of different mixtures 
of four analytes (hydroquinone, tryptophan, phenylalanine 
and DOPA) that were measured in a Perkin-Elmer LS50 B 
fluorescence spectrometer. The measurements of emission 
spectra at multiple excitation wavelengths give rise to an 
excitation-emission matrix (EEM) for each sample and thus 
form a degree 3 tensor which is known to obey the trilinear 
model where the rank is determined by the number of fluo- 
rophores ifTTl [T8l ITT1 . As described in |19| this data set is 
highly suited to assess the performance of different meth- 
ods due to its realistic noise from the physical environment 
and the availability of a priori knowledge of the underlying 
components. Figure [2] shows the estimated emission and 
excitation spectra for PARAFAC with 4 and 5 components 
and AROFAC2 with auto-detected 5 components. We note 
that the results of AROFAC2 are in excellent agreement with 
the known spectra (cf. Figs. 2 and 3 in ifTTlD while the com- 
ponents found by PARAFAC lack accuracy. In addition we 
observe that the AROFAC2 loadings better fulfill the non- 
negativity constraints even though they were not enforced 
explicitly. Furthermore, the peak of the fifth component 
around 315 nm in both excitation and emission spectra which 
can be attributed to Rayleigh scatter in all samples [20; [TT] is 
sharper and thus in better agreement with its expected shape 
when identified by AROFAC2. 

6. CONCLUSION 

With AROFAC2, we have presented an algorithm which can 
determine the CP-decomposition and the rank of a potentially 
noisy degree three tensor. We argue that due to how the 



Emission loadings Excitation loadings 




'250 300 350 400 450 " 240 260 280 300 

wavelength (nm) wavelength (nm) 



Fig. 2. Emission (left) and excitation (right) spectra for the 
Dorrit data using PARAFAC with r = 4 (top row) PARAFAC 
with r = 5 (middle row) and AROFAC2 with r = 5 automat- 
ically detected (bottom row). 

algorithm is constructed, spurious and unstable components 
in the decomposition are not found; since the convergence 
criterion intrinsically enforces stability and thus informative- 
ness of any found component. Our simulations demonstrate 
that AROFAC2 is competitive to the state-of-the-art method 
PARAFAC, but without the need to provide the true rank 
in advance. Also, AROFAC2 outperforms PARAFAC on 
the Dorrit data set which shows that AROFAC2 is a method 
which is more stable to outliers and the influence of non- 
Gaussian noise. 

AROFAC2 uses the intrinsic algebraic structure of a low- 
rank degree tensor in the calculations, as opposed to most 
standard methods such as PARAFAC which assume a model 
and try to fit it, agnostic of its inner structure. We thus argue 
that algorithms exploiting this structure are to prefer when- 
ever available, and the proper starting point for any method 
approaching any problem with algebraic features. We em- 
phasize the potential benefit from applying structural insights 
to construct structure-aware methods. 
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